#### Simulation AR(1) and AR(2) for testing the first order autocorrelation
#### In sample properties
#### Generate Table 1 and 2
#### The user can change vecrho for the autocorrelation values and vecrho1 for the alternatives used for the roots
#### alpha, the significance level, vecn, the vector of sample sizes, nsimu the number of simulations drawn are also global parameters fixed by the user

### Library #####
library(MASS)
###Rho
vecrho=c(0.2,0.5,0.7,0.95,0.97,0.99)

### Output
nom=paste("Table1-2.out")
write("ALL the results in one file",file=nom,append=FALSE)

## Significance level
alpha=0.05

##Sample sizes
vecn=c(50,100,250,500)
sizen=length(vecn)

## Number of replications
nsimu=10000



##### Function to calculate the different p values of the different tests
testcor<-function(n,nns,rho11,rho12){
  buffer=100
  #n:size
  set.seed(nns)
  eps=rnorm(n+buffer)
  denvar= ((1-rho12)^2-rho11^2)*(1+rho12)/(1-rho12)
  y0=rnorm(1)/sqrt(denvar)
  ym1=rnorm(1)/sqrt(denvar)
  
  ytsim=filter(eps,c(rho11,rho12),method="recursive",sides=1,init=c(y0,ym1))
  yt=ytsim[(buffer+1):(n+buffer)]
  eps=eps[(buffer+1):(n+buffer)]
  ytm1=ytsim[(buffer):(n+buffer-1)]
  
  hatrho=mean(ytm1*yt[1:n])/mean(ytm1^2)
  epst=yt[1:n]-ytm1*hatrho
  mt=(epst[2:n]*epst[1:(n-1)]-(1-hatrho^2)*epst[2:n]*ytm1[2:n])/(sqrt(hatrho^2)*mean(epst^2))# our moment on estimated residuals
  mtst=mt/sqrt(var(mt))
  
  ############ Correlation ############
  empcor<-function(hh){
    epstmh=epst[1:(n-hh)]
    epstloc=epst[(hh+1):n]
    return(cor(epstmh,epstloc))
  }
  ######################################
  
  #########
  bp<-function(h){
    veccor=sapply(seq(1,h,1),empcor)
    #### Box Pierce ######
    bp=n*sum(veccor^2)
    rejetbp=(pchisq(bp,h-1, ncp = 0, lower.tail = FALSE, log.p = FALSE)<alpha)
    ####Llung Box #####
    poids=(n-seq(1,h,1))
    lb=n*(n+2)*sum((veccor^2)/poids)
    rejetlb=(pchisq(lb,h-1, ncp = 0, lower.tail = FALSE, log.p = FALSE)<alpha)
    
    return(c(rejetbp*1,rejetlb*1))
  }
  ####My test #####
  rejetmt=(pchisq(n*mean(mt)^2,1, ncp = 0, lower.tail = FALSE, log.p = FALSE)<alpha)
  rejetmtst=(pchisq(n*mean(mtst)^2,1, ncp = 0, lower.tail = FALSE, log.p = FALSE)<alpha)
  matoutcome=c(bp(2),bp(3),rejetmt*1,rejetmtst*1)
  return(matoutcome)
}
####################################################
# Proc used to calculate the rejection frequencies
rejfreq<-function(rho1,rho0){
  rho11=rho0+rho1
  rho12=-rho0*rho1
  
  
  ####### We do it in parallel ####################################
  library(foreach);
  library(iterators);
  library(doParallel)
  nworkers <- detectCores()
  cl <- makePSOCKcluster(nworkers-1)
  clusterSetRNGStream(cl, c(1,2,3,4,5,6,7))
  registerDoParallel(cl)
  


  matsortie1=matrix(NA,6,sizen)
  for (ln in seq(1,sizen,1)){
    matrej<- foreach(ns=1:nsimu,.packages = "MASS", .export=c("testcor","alpha","vecn","sizen"), .combine=rbind) %dopar% {testcor(vecn[ln],ns,rho11,rho12)}
    matsortie1[,ln]=apply(matrej,2,mean)*100
  }
  stopCluster(cl)
library(xtable)
write(paste("RESULT DIFFERENT TESTS ","- rho0=",rho0,"- rh01=",rho1),file=nom,append=TRUE)
rownames(matsortie1)=c("BoxP-2","LlungB-2","BoxP-3","LlungB-3","Robust","RobustIS");
colnames(matsortie1)=vecn
print(xtable(matsortie1,caption=paste("Autocortest-","rho0=",rho0,"-rho1=",rho1)),digits=2,file=nom, append= TRUE)
return()
}  
#####################################################################
  
  

#### One big prog to calculate the same quantities for different values of rho0
  
bigloop<-function(rho){
#### Size and power through the choice of rho1
vecrho1=c(0,0.1,0.2,0.3,0.5)
sapply(vecrho1,rejfreq,rho0=rho)
}


### Table 1 and 2 outcomes + more
sapply(vecrho,bigloop)
